Parametric resonances in electrostatically interacting carbon nanotube arrays 
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We study, numerically and analytically, a model of a one-dimensional array of carbon nanotube 
resonators in a two-terminal configuration. The system is brought into resonance upon application 
of an AC-signal superimposed on a DC-bias voltage. When the tubes in the array are close to each 
other, electrostatic interactions between tubes become important for the array dynamics. We show 
that both transverse and longitudinal parametric resonances can be excited in addition to primary 
resonances. The intertube electrostatic interactions couple modes in orthogonal directions and affect 
the mode stability. 



I. INTRODUCTION 

During recent years, several experimental realizations 
of nano-electromechanical (NEM) resonators based on 
carbon nanotubes (CNT) or carbon nanofibers (CNF) 
have been demonstratedii^i^ii^i^i^. Carbon nanotubes 
have established themselves as strong material candi- 
dates for use in NEM-resonator systems, partly due to 
their favorable mechanical properties^, such as low mass 
and high elastic modulus. Thus, using CNTs/CNFs al- 
lows for operational frequencies of NEM-resonators that 
reach into the GHz regime. In the most recent experi- 
ments, the long predicted high quality factors of the order 
of Q ~ 10'^ have finally been achieved. This makes these 
resonators interesting from a technological point of view 
with application areas such as tunable RF-filters and fast 
low power switching elements^'^'^. However, for such ap- 
plications, a major drawback is the high impedance lev- 
els offered by single nanotube devices. The ensuing low 
power transduction^ makes integration of such devices 
with current state-of-the-art CMOS technology difficult. 
One way to overcome this problem is to construct de- 
vices based on parallel arrays. For arrays it is desirable 
to know how interactions betwen elements affect the op- 
eration of devices. It is for instance important to know 
how closely spaced array members can be placed without 
drastic changes in performance. 

Apart from the technological incentive to study 
MEM/NEM arrays, the problem is also of fundamental 
interest. The dynamic response of coupled NEM/MEM- 
resonator systems in combination with nonlinearities is 
known to lead to unexpected and/or unintuitive behav- 
ior. Examples are intrinsic localized mode a^^i^^'^'^i^^ and 
mode synchronizations^'^®. In addition, parametric reso- 
nances in MEM/NEM arrays have been studied both ex- 
perimentally and theoretically-'iii^'iSii^SiSi . In those stud- 
ies parametric response was induced by applying an AC- 
voltage component between alternating beams in beam- 
arrays. Parametric resonances can be narrower than fun- 
damental resonances and finds use in for instance para- 
metric amplifiers. 

In this paper we study theoretically a vertical one- 
dimcnsional regular array of CNT/CNF resonators (see 
figure 1). When the resonators are not too widely 
separated, electrostatic interactions between the tubes 




Figure 1: (Color online) Schematic layout of a vertically 
oriented carbon nanotube NEM array. Through patterning 
of catalysts, regular one dimensional arrays of carbon nan- 
otubes or nanofibers may be grown on top of an electrode 
(Source). A second electrode (Drain) is formed from a de- 
posited layer at the same height as the tube tips. Actuation 
and transduction of tube displacement is achieved through 
electrostatic (capacitive) coupling between the tubes and the 
drain electrode. 



become important and affect the dynamical response 
of the system. While the system considered in this 
paper shares some features with previously studied 
systemsiii2iii^i22ii^, there are several important differ- 
ences. Among them two are worthy of special attention. 
Firstly, all the tubes in the system are connected to the 
same voltage source, containing both a DC- and an AC- 
component, thereby eliminating the need for individual 
contacting of alternating array members. This makes 
the interactions between tubes repulsive, rather than at- 
tractive. Secondly, in contrast to beams with rectan- 
gular cross-sections, where the characteristic vibration 
frequencies differ for vibrations in different directions, 
CNT/CNF have a circular cross-sections and motion in 
two dimensions plays an important role. 

We find that the fundamental resonance, where the 
tubes oscillate in unison towards the drain electrode, is 
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not drastically affected by interactions. However, several 
new resonances appear including longitudinal resonances, 
where the tubes oscillate in the direction along the ar- 
ray. For small arrays, these resonances have the form 
of hardening Duffing type resonances which develop into 
a band of resonances as the array gets larger. Further, 
two parametric resonances are present in the systems. 
Both transverse as well as longitudinal motion may be 
parametrically excited, both with multiple branches and 
complex bifurcation structures for the larger arrays. The 
electrostatic coupling between tubes also affect the sta- 
bility of the longitudinal motion which becomes unstable 
due to parametric excitation of transverse oscillations. 

We begin this paper by presenting, in section [ill a- 
simplified lumped electromechanical model to derive the 
main qualitative features of the system. Then, in sec- 
tion mil based upon numerical integration of the equa- 
tions of motion, the general features of the dynamic re- 
sponse of a one-dimensional CNT resonator array are dis- 
cussed. To better understand the main characteristics of 
this response we focus on the the smallest possible array 
(two tubes) in section HVl The two-tube system is treated 
both numerically and analytically. We derive frequency 
response equations and analyze stability for the various 
resonances using perturbation theory. These analytical 
results are found to work as good approximations for find- 
ing the loci of the resonances also in larger systems. This 
treatment is then followed up, in section |Vl with a dis- 
cussion of how the response changes as the arrays become 
larger before concluding in section IVll 

II. MODEL 

For an array consisting of N tubes we denote the co- 
ordinates of the central axis of each (undeformed) tube 
by X" = {Xi,Yi) where i = 1,...,N. For the vibra- 
tions of the beams we consider only excitations of the 
fundamental flexural modes for which a lumped model 
is suitable^i22,. Describing the position of the tip of the 
cantilevers by the coordinates we use the equations 
of motion—: 

m,X, + ma±; + m^cul,{X, - X^) = Pf ^ = 1, .., N. 

Here are the effective masses of the tubes and cuQi 
the natural resonance frequencies. For tubes with cir- 
cular cross sections these values ar e given b y^^ nii — 
pAiLi/5.68i and ujoi = 3. 516 ^ ^/Eli/piAi where Ai 
is the cross-sectional area, li the moment of inertia, E 
the Young modulus and Li the length of tube number 
i. We have also introduced a viscoelastic damping term 
7i for each tube to account for mechanical losses. In 
the absence of a gaseous medium surronding the tubes, 
the main sources of dissipation are clamping losses and 
Ohmic losses. 

Aside from elastic forces, external electrostatic forces 
Ff- act on the tubes. These forces depend on the geom- 
etry as well as the instantaneous charge distributions on 



the tubes. To find the exact charge distributions on the 
tubes is a hard problem. Numerical simulations using 
FEM and BEM of electrostatically interacting tubes^— 
reveal that the charge is mainly located at the tip of the 
tubes. Furthermore, the main contributions to the bend- 
ing moments arise from forces close to the tube tips. Thus 
to model the electrostatic forces we make the simplified 
assumption that the charge Qi on each tube is concen- 
trated to a conducting spherical shell located at the tip 
of each tube. The model we consider is one of conducting 
spheres, attached with springs to their equilibrium posi- 
tions and being able to move in the a;y-plane. The drain 
electrode, which can be taken to be at zero potential, is 
modeled as an infinite conducting plane. There is also 
a possibility of actually having metallic grains on top 
the tubes. In plasma-CVD growth of CNT/CNF from 
Ni catalysts, tip-growth results in the metallic catalysts 
residing at the tips of the tubes. 

In general, the full charge distribution, rather than the 
total charge, on each sphere is needed to find the forces. 
Provided the intertube separation as well as the tube- 
drain separations are larger than the tube diameters, the 
dominant contribution to the electrostatic forces comes 
from the monopole contributions of these charge distribu- 
tions. Restricting attention to this case the electrostatic 
free energy of the system [the tubes being biased by the 
common time-dependent voltage V{t)] is 

N 

J'cl. — J'scU + J'int. — V{t) ^ Qi, 

i=l 

where the electrostatic self-energy is 

^ i=i * 

and the total electrostatic interaction energy is 

.7-1 \ ^ \ ^ QiQj 1 \ ^ QiQj 

~ 4^ ^ ^ |X, - X,r 8^ |X^ - X,| ■ 

Here X^ = {—Xi, Yi) denote the images of X^ = {Xi, Yi) 
in the plane X — 0. The charge distribution in the array 
is then found from solving the linear system dJ-'d. / dQi — 
after which the electrostatic forces may be found as 
Fi = — Vxi.?^ci.- 

We consider now a uniform system with identical 
tubes, i.e., Di = D, ujoi = ujq, 7i = 7 and = m. 
Rescaling the coordinates to dimensionless form accord- 
ing to Xi = Dxi, Yi = Dyi yields 

x,+7X,+a;2(^_^0)^J_Fei._ 

In the same way we rescale electric quantities through 
introducing a unit voltage Vq, i.e. V = vVq, Qi — 
qiAneoDVo and the electrostatic charging energy £c = 
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Figure 2: (Color online) Time averaged mechanical energy 
for a harmonically driven array with 100 tubes. The blue 
solid lines represent the mechanical energy stored in trans- 
verse modes (vibrations towards the drain electrode). The 
red, dashed lines represent the energy stored in longitudinal 
modes (vibration along the array axis). Arrows indicate the 
directions of transitions in the hysteretic region. Top: Inter- 
tube separation 200 nm. Middle: Intertube separation 150 
nm. bottom: Intertube separation 125 nm. 



AttcoDYq. The corresponding dimensionless electrostatic 
free energy /d. = Td./Sc is then 



III. RESPONSE TO HARMONIC DRIVING 

Figure [5] shows the response of an array with 100 
tubes when it is driven with an AC-signal on the source 
in combination with a static DC-bias voltage V{t) = 
Vo + Vi cosiuDt) (DC-bias voltage Vq = 12 V, AC-signal 
Vi — 2 V). The tubes are each 1 fim long with a diam- 
eter of 25 nm. The distance to the drain electrode is 
150 nm. Plotted are the two orthogonal components of 
the mechanical energy corresponding to transverse (blue 
solid line) and longitudinal vibrations (red dashed lines). 
The energy is scaled in terms of the dimensionless units 
introduced above with a timescale set to to = 0-1 ns. 

In the top panel the spacing between the tubes is 200 
nm. Clearly visible is the primary transverse resonance 
(blue lines) where all the tubes oscillate in phase with 
each other. This mode corresponds to the resonance of 
a single nanotube. A band of longitudinal modes can be 
seen (dashed line) just above 100 MHz. In the middle 
panel the tubes are more closely spaced (150 nm) and 
two additional resonances are present. The transverse 
(around 190 MHz) is a parametrically excited resonance 
where each tube oscillate with half the driving frequency. 
In this resonance neighboring tubes oscillate with oppos- 
ing phases (optical mode). Above 210 MHz, is another 
parametric resonance in the more closely spaced arrays. 
This is a parametric resonance of the band of longitudi- 
nal modes of the array. In the bottom panel the spac- 
ing has been narrowed down further to 125 nm. The 
parametric resonances are now stronger and the primary 
transverse resonance has become hysteretic. The appear- 
ance of hysteresis can here be understood by considering 
the attractive force between a tube and the drain elec- 
trode. For widely separated tubes, each tube is attracted 
by its own image potential alone, while for more closely 
spaced tubes, the images charges from neighboring tubes 
contribute to this force. 

In the next section we show how these resonances and 
their main characteristics can be understood from an- 
alyzing a two-oscillator array. Then, in section |V] we 
study, numerically, how the response changes qualita- 
tively as the size of the arrays grow larger. The para- 
metric resonances, both the transverse and the longitu- 
dinal, show a complex behavior with multiple bifurcation 
points. These are not shown in the panels of figure[2]but 
will be addressed further in section IVl 
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In order to study the response of the system we solve 
the dynamic equations ([1]) numerically using a velocity- 
Verlet algorithm. In the next section we present the main 
qualitative feutures of the mechanical response of the sys- 
tem to a harmonic driving field. A more quantitative 
discussion is then carried out in sections [TVl and fVl 



IV. TWO OSCILLATORS, CASE STUDY 

In this section we study the simplest case, namely an 
array consisting of only two cantilevers with identical 
physical parameters. This case can be analyzed ana- 
lytically and serves to validate numerical modeling and 
provides insights for larger arrays. 

We take the positions of the undeflected tubes to be 



(a;o,yo), 



{xo,—yo) where Xq < and i/q > 



and the drain to be the plane a; = 0. Introducing the 



4 



variables x± = xi ± X2 and y = yi — 1/2 we have the 
equations of motion 



x+ + 7i+ + i^o{x+ — 2xq) = i^ujqv'^ 



9l 



9l 



{x+ — x^Y 



- 251 52 



{x+ + x^)'- 



(x^+y2)3/2 



25152 



9l 



9l 



[x^ + 



y + 72/ + Wo(y + 2?/o) = 2:^a;oW^5i52 



+ y2)3/2 

(xi +y2)3/2 

2/2)3/2 



(2) 



(3) 



(4) 



The relevant electromechanical coupling constant is = 
^c/'^-mcch. = Ec/mD^LOQ. In terms of numbers v « 
lO^'^V^^L^-D"^ if length is measured in nm and the tubes 
are assumed solid with a Young modulus of 1 TPa3&. 
The functions 51,2 are found by solving exactly the elec- 
trostatic problem and are given by 

51 = (.T++x_) + '-^ 



52 



{x+ - x_) 



-x2_)(A2 -4) -4a:+ - 1^ 
(a;++a;_)(A-2)-l 
- .x2_)(A2 -4) -4x+ - 1 



where A is defined as 
1 



x'^_ + 2/2 



+ 2/2 



A_ - A, 



A spectrum that reveals the most important features of 
the response to a harmonic AC-drive on the gate is shown 
in figure [3l This figure was obtained from numerical in- 
tegration of the dynamic equations for a system with the 
following parameters: tube diameter D — 2b nm; tube 
lengths L = 1/im; bare quality factor kq = 100; tube 
positions (Xo,lo) = (—150,62.5) nm; Young modulus 
E = \ TPa; tube density; p = 1.2 g/cm^. The apphed 
voltage to the system was Vb = 12 V and Vi = 2 V. 

The figure was obtained by sweeeping the drive fre- 
quency both upwards and downwards. On the vertical 
axis of figure [3] the dimcnsionless average mechanical en- 
ergy of the tubes is shown. The motions in the longitu- 
dinal direction (2/-direction) and the transverse direction 
(x-direction) have been separated for clarity. Both the 
transverse response (blue line) as well as the longitudi- 
nal response (red) show three main peaks each. We have 
labeled these peaks A,A',C and B,B',D respectively. 
Hysteresis in the frequency plane is present in the peaks 
B, C and D (for the peak C the hysteresis is too nar- 
row to be clearly seen in figure [3]). In the subsections 
below we treat each of these resonances in more detail. 
Note that the subsection labels follow the labelling of the 
peaks in figure [31 
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Figure 3: (Color online) Average mechanical energy of a two- 
tube system in response to harmonic driving. Blue solid 
lines: Average mechanical energy stored in vibrations in the 
x-direction. Red dashed lines: Average mechanical energy 
stored in vibrations in the {/-direction. Both upwards and 
downwards sweeps in frequency were made. A, A': Funda- 
mental resonances in the x-direction. Both tubes oscillate in 
phase towards the drain electrode (a;+-resonance). B, B': 
Fundamental resonances in the {/-direction. The tubes oscil- 
late with opposing phases in the direction parallel to the drain 
electrode ({/-resonance). C: Parametric resonance in the x- 
direction. The tubes oscillate with opposite phases at half 
the driving frequency in the direction towards drain electrode 
(a;_-resonance). D: Parametric resonances in the y-direction. 
The tubes oscillate with opposing phases at half the driving 
frequency in the direction parallel to the drain electrode (j/- 
resonance) . 



We begin by determining the stationary points. For 
small deflections around equilibrium it is sufficient to 
keep only the dominant terms in I/xj^. and 1/y which 
yield the new dynamic equations: 



x+ + 7i+ + ujI{x+ — 2xo) 



y + jy + ujl{y + 2yo) = v- 



,2„,2 



,2„,2 r 



y2)3/2 



(4 



2/2)3/2 



(5) 



In the limit of large intertube separation [y — > 00) one 
retains the result of noninteracting tubes whereas the 
limit X 00 reduces the problem to one in the y- 
direction only. The equations also decouple in the limit 
of small |a;+| (recall that both x+ and y are negative) 
due to screening of the electrostatic interaction between 
the tubes by the drain electrode. The system (O can be 
used to determine the stationary deflections x^ (t) = Xg , 
y(t) = ys in the absence of an AC-component. These 
time-independent solutions are found by solving the sys- 
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tern 



Xs = 2x0 + - 



x2 2/2)3/2 



~2yo + ^ 



.y2)3/2 



(6) 



For small biases vq (i.e. <C 1) we solve pertubatively in 
v and get 



2xa 



1 



,-3 



1 



16yg 

V 



3 V'O 



So) 



(1 - So) 



(7) 



where co = \xqIvq\ and Eo = (1 + o^o)^^^^- The term 
So accounts for mutual screening of the tubes (if do — 
l^o/yo| ~* then So — > 1). This approximation is valid 
for small static deflections and far away from snap-in (the 
tubes making contact with the drain). Note that the 
system ([71) rests on the approximation |l/a;+| <C 1 and 
1 1/2/ 1 ^ 1- This is consistent with the approximation 
of only keeping the monopole contribution to the total 
charge distribution on the tube tips. 



A. Fundamental transverse resonances, 
(a;+-resonance) 

In the fundamental transverse resonances [A and A! in 
figurc[21) the tubes oscillate in phase with each other. Nei- 
ther of these resonances differ appreciably in nature from 
those of single tube systems. The subharmonic arises 
from the double frequency component of the driving oc- 
curing due to the term. The main resonance tunes 
downwards in frequency with increasing bias and has a 
Dufflng type nonlinearity of the softening kind. We omit 
the analysis of this resonance in this paper since it has 
been already thoroughly studied previously in the littera- 
ture in conjunction with single cantilever resonators (see 
for instance Ref. [24 



B. 



Fundamental longitudinal resonance, 
(j/- resonance) 



The fundamental y-resonance (resonance B in figure[31) 
has the shape of a Duffing resonance with a hardening 
nonlinearity. A more clear view of the resonance is seen 
in the inset of flgure[31 In this figure a false color plot of 
the mechanical energy in the primary longitudinal reso- 
nance (y-resonance) as a function of bias voltage (vertical 
axis) and drive frequency (horizontal axis) is shown. The 
figure was created sweeping the frequency downwards. 
Clearly visible is the upwards tuning of resonance fre- 
quency with increasing bias and the sharp onset of res- 
onance at the bifurcation point. The inset shows the 
response along the 9V bias cross-section (along the dot- 
ted line). The thick black curves are the results from 
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Figure 4: (Color onfine) False cofor plot of mechanical en- 
ergy of the primary fongitudinal resonance (j/-resonance) as 
a function of bias voltage (verticaf axis) and drive frequency 
(horizontal axis). The figure was created sweeping the fre- 
quency downwards. Tlie solid white line corresponds to 



the expression (|12p for the central frequency 



Sweeping 



frequency downwards results in an abrupt change in the re- 
sponse at the bifurcation frequency oj^' (the high-frequency 
edge of the large ampiitude (red) region) . The dashed white 
line shows the the bifurcation point found from solving the 
frequency response equation (I13|) . The inset shows the re- 
sponse afong the cross section at a bias of 9 V (aiong the 
dotted line). The thick bfack curves come from numericai 
simulations whereas the dash-dotted blue lines are soiutions 
to the frequency response equation (|13|l . 



numerical simulations of the dynamic equations ([T]) . The 
solid black line correspond to downard frequency sweep 
while the dashed line to upward frequency sweep. 

The characteristics of the fundamental y-resonance can 
be found using perturbation theory. Considering this res- 
onance we take a;+ = Xg and assume x_ = 0. When 
x_ = the product (71(72 simplifies to 
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[A. 



„-ll2 ■ 



(8) 



For small oscillations we may further set 



■ V 



Vs 



a. 



(9) 



To obtain an estimate of the parametric dependence 
of the resonance character we analyze the system us- 
ing the method of averaging^^ by making the Ansatz 
y{t) = ys + Y{t) cos{uiDt) in response to a drive given by 



[1 -|- 2ecos{uiDt + S)]. Assuming Y/Y ujo, 
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Figure 6: (Color online) Fundamental y-resonance along the 
dashed line of figure [l] Blue solid lines were obtained from 
solving equation l|13p while the black dashed lines are from 
immerical simulations. Above the red dash-dotted line, per- 
turbation theory predicts the fundamental longitudinal reso- 
nance to be unstable towards parametric excitation of the x-- 
resonance. The point of destabilization of the upper branch 
occuring where the dash-dotted red line and the upper blue 
line cross. 



Figure 5: Mechanism for destabilization of the high ampli- 
tude branch of the fundamental y-resonance (B-resonance in 
figure |3]). The figures show the time evolution of the trajec- 
tories of the tubes in the a;j/-plane during destabilization. A: 
Initially the tubes are executing high amplitude motion in the 
y-direction. The coupling to the motion in the x— direction 
causes instability of the X- mode which starts to grow (B). 
After having reached high amplitude in both x— and y— di- 
rections (C) irregular motion ensues (D) before the motion 
finally settles in the low amplitude branch of the y-resonance 
(F). 



the differential equation for the amphtude Y is 



Vs 



Y = ujoiy^sm6[K2{Y,ys)-Ko{Y,ys)]-^jY 



Yujd = Yuq 
Here 

Kn{Y,ys) = 



1 



uJoiy{Ki+ecos5[Ko+K2]). (10) 



d(j) cos ncf)— 
TT ./o (y. 



FC0S( 



Ay 



(11) 



with A EE + 

For small oscillation amplitudes where the response 
does not bifurcate we solve the system to first order in v 
in the limit Y ^ 0. For the amplitude Y and the center 



of resonance Wc*''' we get 



Y 



ujQiye (1 - Eo) 



,{y) 



Jy^/ujo = 1 + 



16%' 



1 



2^0 -lis, 



^5/3 



(12) 



Equation is useful for estimating the center fre- 

quency even for large oscillations as can be seen in fig- 
ure m where uj^^ obtained from equation p2p is drawn 
as the solid white line. As the drive gets stronger bifurca- 
tion occurs (see figure |3]) and the shape of the resonance 
is found from solving the frequency response equation 



2 2 2 



^2y2 



[K2 - K^? 



[2Y{lod - cjq) - ujqvKi 



(13) 



The blue dash-dotted lines of figure 2] depict the reso- 
nances obtained from equation (jl3[) , showing good agree- 
ment for small amplitudes. The locus of the bifurcation 
point ujg may also be determined from solving equa- 
tion (fT3|) and this solution for w^'' is shown as the white 
dashed line in figure U) 

Whereas both the location of the resonance and the 
bifurcation point can be estimated using Eq. (fTU|) . this 
is not true for finding the extent of the hysteresis. The 
destabilization of the high amplitude branch is connected 
with an instability towards resonance of the 2;_-mode. 
The process of destabilization is depicted in figure [H] 
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where the entire time evolution of the trajectories of the 
tubes are shown. Starting in panel A, the system is in 
the high amplitude branch. The instability towards res- 
onance of the a;_-mode causes an increase in motion in 
the transverse direction (panels B and C). Decay to the 
lower branch (panel F) occurs through irregular motion 
of the tubes (panels D and E). The location where the 
high-amplitude branch of the y-resonance becomes unsta- 
ble can be found analytically using perturbation theory. 
In figure [S] the region of instability towards parametric 
excitation of the a;_ -resonance is shown as the red dash- 
dotted line. As can be seen, a good estimate of the locus 
of the destabilization can be determined. The perturba- 
tive analysis is found in Appendix O 



C. Instability towards parametric resonance 
(a;_-resonance) 

We now turn the attention to the parametric resonance 
of the a;_-mode (the C-resonance of figure [3]). Writing 
out explicitly the right hand side of the equation of mo- 
tion ^ we have 

x^+^x^+Lolx^ ^ ~2x^LUQ[l + 2eco8{LLiDt)]F{x^) (14) 

with 

2x^(A-2)2-2(A-2) 



F = 1^ 



5152 



[(x2 - a;^)(A2 - 4) - 4xs - l]^ {x'i_ + yl)^!"' 



The right hand side is proportional to x- characteristic 
for a parametric drive. A simple parametrically driven 
harmonic oscillator 

X + "fx + = ^^(^^[l + 2e cos(cij£)t)]a; 

will be unstable^'' if uj'^ < ujd < i^t, where 



UJQ 



if2 



the 



provided 



l-K 

(15) 

discriminant is positive, i.e. 
To find the point of instabil- 



ity we keep only the lowest order term in a;_ in 
equation (fT^. This gives 



K = 2u 



(2-Ao)(2+_a3)x, + al 
2;2(Ao + 2- 



-1^2 



(16) 



where Xs and ys are the stationary points and as = Xs/ys 
and Aq = ^A+ — y^^. A comparison between numeri- 
cal simulations and the region of instability is shown in 
figure [T] For small biases the agreement between theory 
and numerics is good while it deviates for larger biases. 
This deviation is due to the approximate relations ([7]) to 
find Xs and ys- 

For larger amplitudes we must consider the full equa- 
tion of motion p4|) . Introducing action angle coordinates 
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Figure 7: (Color online) False color plot of the mechanical 
energy in the parametric a; --resonance (resonance C in fig- 
ure |3]) as a function of bias voltage (vertical axis) and drive 
frequency (horizontal axis). The dashed white line denotes 
the region of instability according to equations ([T3|l and fTO]) . 



X- = X {t) cos[(t){t)] and i_ — — X(i)w sin[0(i)], and ex- 
panding _F in a Fourier series 

F{x-) = F{X COS0) — a(j/2 + ^ a„ cosncj), 

we obtain after averaging out fast variables the au- 
tonomous system 



X ^ X 



wg(a4-ao) . 1 

e— smd 7 

2u} 2 ' 



ujQ{ao + 2a2 + 04) 



q) — Lo + e 
u} — WqVI + ao + a2 



2u> 



cos (5 



Here 5 is the relative phase of oscillation with respect 
to the drive. We note that in the limit X — > we have 
0-2, ^4 = and gq — K. A comparison between the results 
of perturbation theory and numerical simulation is shown 
in figurelH Here the mechanical energy in the parametric 
^--resonance is shown for a bias of Vq = lAV (black solid 
line is the downward frequency sweep and dashed line the 
upward frequency sweep). The red dash dotted line is the 
result of solving the frequency response equation 



,2^0 



7 



(tOD - 2ujf 



uP- (04 — ao)2 (ao -f 2ai + a4)2 



(17) 



While agreement between perturbation theory and nu- 
merics is good it does not work well close to the point 
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Figure 8: (Color online) Close-up of the response of the 
parametric a;_-resonance (Cross-section at Vb = 14 V of fig- 
ure [Tj. The resonance has the characteristics of a parametri- 
cally driven DufBng resonator with a hardening nonlinearity. 
Near the point of instability of the upper branch, coupling 
to the longitudinal mode causes beats where energy is trans- 
ferred between transverse and longitudinal modes periodically 
in time. Black curves are from numerical simulations while 
the red dash-dotted curves come from solving the frequency 
response equation (|17p . The inset shows the average me- 
chanical energy stored in transverse and longitudinal modes 
as function of time. The red curve, showing energy for the 
longitudinal mode, has been magnified 500 times and verti- 
cally displaced for clarity. 



of instability of the upper branch. Here, there is noise 
in the curve obtained from numerical simulations. This 
noise comes from coupling to the longitudinal mode. The 
inset in figure [8] shows how the average energy stored in 
transverse (blue) and longitudinal modes (red) vary in 
time. The curves have been displaced for clarity and the 
red curve is magnified 500 times. While the energy trans- 
ferred to the longitudinal mode is very small compared 
to the energy in the transverse mode, the excited lon- 
gitudinal vibrations has great impact on the transverse 
vibrations. 



D. Parametric longitudinal resonance (j/-resonance) 

Finally we study the conditions for observing the para- 
metric longitudinal resonance (D- resonance in figure [31) ■ 
As in the case of the parametric resonance in the trans- 
verse direction, the region of instability in the frequency 
plane towards parametric resonance in the y-direction is 
determined by the equation (jlSp. Starting from the equa- 
tion ([4]) and making again the approximations in ([8|) and 
we find 



K 



{ys ~ Af 



V 



1 + -s. 

2 



(18) 



Figure 9: (Color online) False color plot of the mechanical 
energy in the parametric longitudinal resonance (y-resonance) 
as a function of bias voltage (vertical axis) and drive frequency 
(horizontal ajcis). The dashed white line denotes the region of 
instability according to equations (|15|l and (|18|l . The inset 
shows the response along the cross section at 9 V bias (dotted 
line) . The thick black curves come from numerical simulations 
and the solid blue lines are solutions to the frequency response 
equation (|19p . 



A comparison between numerical simulations and the re- 
gion of instability is shown in figure [H The figure was 
created sweeping the frequency downwards and the bifur- 
cation edge is visible as the sharp transition between dark 
(blue) and bright (red). For small biases the agreement 
between theory and numerics is good while it deviates 
for larger biases. This deviation is again due to using the 
approximate relations ([7]) to find Xs and ys respectively. 

As in the preceeding subsections we may use per- 
turbation theory to study the large amplitude response 
of the parametric resonance. Assuming yit) = ys + 
Y coslujBt/2) and v{t)'^ « v'^{\ + 2ecos[ujDt -\- 5]) the fre- 
quency response equation can be derived 



yMi + ^-lJ^' 



(Ki + Kz) 



tan 5 



(19) 



where Kn are given by equation (jlip . A comparison 
between perturbation theory and numerical simulations 
is shown in the inset of figure [H Again agreement is 
good but fails to predict where the upper branch be- 
comes unstable. The destabilization of the parametric 
y— resonance occurs in the same way as the fundamental 
y— resonance, i.e. through parametric excitation of the 
2;_-mode and can be analyzed following the along the 
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Figure 10: (Color online) Response of a four-tube system 
with the same physical parameters as the one in figure [S] The 
inset shows a close-up of the parametric transverse resonance 
with the directions of the transitions in the frequency plane 
indicated by arrows. 
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Figure 11: (Color online) Response of an 8-tube system with 
the same physical parameters as the one in figure [3] The in- 
set shows a close-up of the parametric transverse resonance 
with the directions of the transitions in the frequency plane 
indicated by arrows. In the area denoted noisy region, the 
coupling to longitudinal motion causes the amplitude of trans- 
verse motion to oscillate in time. 



lines of the calculation in Appendix \^ 



V. SEVERAL OSCILLATORS 

Having treated the two-oscillator system in some detail 
we now move on to describe how the system response 
changes with increasing system size. For this we use the 
same system parameters (geometry and bias voltages) as 
those used to obtain figure[3]and only change the number 
of tubes in the array. We have done detailed simulations 
for systems with 4, 8, and 16 tubes and the corresponding 
frequency responses are shown in figures [1011121 
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Figure 12: (Color online) Response of a 16-tube system with 
the same physical parameters as the one in figure[3] The black 
dots indicate that when the parametric longitudinal resonance 
became unstable, strong excitations of transverse modes oc- 
cured that lead to snap-to-contact of the system. The inset 
shows a close up of the parametric transverse resonance. In 
the area around 195 Mhz irregular behavior with high ampli- 
tude motion occurs. 



The fundamental transverse resonance is not markedly 
affected by the increasing array size. This is ex- 
pected since here all tubes oscillate in phase with each 
other. The fundamental longitudinal resonance is how- 
ever strongly affected, the single, hysteretic peak from 
the two-tube system, developing into a broad band of 
excited oscillation modes. The presence of this band is 
reflected also in the longitudinal parametric resonances, 
where the development of band structure is present in 
terms of multiple branches and bifurcations in the re- 
sponse. This type of behavior has been seen in paramet- 
rically driven NEM/MEM arraysi^. Also for the larger 
arrays large amplitude excitations of longitudinal oscil- 
lations can be destabilized due to parametric excitation 
of transverse modes. In figure [12] two particular such 
points are marked with black circles. At these points the 
excitation of the transverse modes became so strong that 
snap-to-contact occured. 

Also the parametric transverse resonance shows the de- 
velopment of a band structure. In contrast to the funda- 
mental resonance where this band structure is not acces- 
sible, several branches can be reached through parametric 
excitation. In figures [T0l[T2l the insets show closeups of 
the parametric transverse resonances. 

While more and more modes appear as the arrays get 
larger, one feature is common to all the systems. This 
feature is the noisy region around 195 MHz. In this re- 
gion, energy is transferred between transverse and longi- 
tudinal modes just as in the case of the two-tube system 
(see figure [8]) but without destabilizing the transverse 
motion. 

As for the location of the resonances in the voltage- 
frequency plane these do not differ appreciably from the 
two-tube system and the perturbative formulas derived 
in the preceeding section can be used to estimate if and 
where the system will be unstable to a certain resonance. 
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VI. CONCLUSIONS 

In order to investigate the effects of electrostatic in- 
teractions between carbon nanotubes in NEM-resonator 
arrays we liave studied a simple model both analyti- 
cally and numerically. We have found that, apart from 
excitation (fundamental and parametric) of a band of 
longitudinal modes, also parametric excitation of trans- 
verse modes is possible. With increasing number of res- 
onators, these resonances become successively more com- 
plicated and exhibit rich behavior with several overlap- 
ping hysteresis loops, bifurcation points etc. The trans- 
verse modes are also responsible for destabilizing the lon- 
gitudinal modes at high amplitudes and may lead to snap 
to contact. Also, the parametrically excited transverse 
modes, show regions of irregular behavior coming from 
coupling between transverse and longitudinal modes. We 
have shown, that the features of the response of ID- 
arrays can be understood qualitatively through study- 
ing the simplest possible array, a two-tube system. Also 
quantitative predictions based on the two tube system 
can be used to obtain estimates of regions of instability 
towards parametric resonances and to estimate frequency 
tuning. 

From a technological point of view, these estimates 
can help in designing array resonator systems to avoid 
unwanted resonances while maintaining a high packing 
density. Utilizing parametric resonances could also be a 
path to further increase the operation frequency in tech- 
nical applications and by tuning the bias voltages the 
width of the region of instability can be tuned to an ar- 
bitrarily narrow frequency domain. So far, only uniform 
arrays have been studied. For applications, disorder must 
be accounted for and further studies are needed. 



Acknowledgments 

This work was supported by the Swedish Foundation 
for Strategic Research (SSF) and the EU through the 
Nano-RF project FP6-2005-028158.This publication re- 
flects the views of the authors and not necessarily those 
of the EC. The EC is not liable for any use that may be 
made of the information contained herein. 



Appendix A: DESTABILIZATION OF PRIMARY 
y-RESONANCE 

We here give a brief derivation of the criteria for desta- 
bilization of the fundamental longitudinal mode through 
parametric excitation of the transverse a;_-mode of the 
two-tube system. Following the same lines, the stabil- 
ity of the parametric longitudinal excitation can be ana- 
lyzed. 



The longitudinal vibrations are destabilized by the x- 
mode, which has the equation of motion 



•• I • I 2 2 2 

X- -\- 72;- -I- UJqX- — VLOqV 



9l 



9l 



{x+ -x_)' 



+ 2.91^2- 



X- 



.y2)3/2 



For small oscillations of the x^-mode, the right hand 
side can be approximated for large amplitudes Y of the 
y-mode (recalling that y — i/s + Y cos ivot) yielding the 
equation 



X- + 7i_ -I- ujqX^ 



-x^F{t) 



where 



F{t) 



2vA^ujIvI{1 + 2ecos{ijjDt + 6)) 
{ys+Ycos{ujDt)){Ycos{LJDt)+ys - A)"' 



(Al) 



and we have defined A = {2 + 1/xs — a)^^ and a = 
{xl+y^)~^^'^ respectively. Changing to action angle vari- 
ables {x- = X{t) cos[(j){t)], X- — —X{t)LJsm[(j){t)]) and 
averaging over fast variables results in the autonomous 
system 



X = X 



(sin 2(l)F{t)) 7 
2 



2^cjg + Co 



UJ0+C0 + 



1 



2y/ujl+Co 



(cos20F(t)) 



where the brackets denotes the averaging (•) = 
(27r)^^ J^^ ■ d(j), and we have expanded F{t) in a Fourier 
series F{t) — c„e™'^°*. At the onset of the destabiliz- 
ing ^--resonance we have 4> = lod and (j) = i^ot+'d- Eval- 
uating the averages and setting C2 = |c2|e*''' one finds: 



X ^ X 



|C2| 



UJD 



2y^ 



I sin(2z9 - A) 7 
2v/^o + Co 2 



:|c2|cos(2i9- A). 



Co 



The region of driving frequencies where the high ampli- 
tude branch of the y-mode can be destabilized by the x_ 
mode can then be found as uj 
with 



dcstab- 
D 



, ^ destab.+ 



UJ 



dcstab.ib 
D 



- Co ± 



C2 



4(t^g + Co) 



Using the expression (|A1[) the Fourier coefficients cq and 
C2 can be evaluated exactly: 
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Co 



= 2v 



Y + 2{A- ys)ecos5 r + 2ey,, cosJ A{A- ys) + 2AYecos5 



Y^-Y^ + {A~y,f Y^-Y^+yl {-Y^ + {A - y^ff'^ 



(A2) 



C2 



2v 



(y2 + _ y^^y^ ^(2F)e cos 5)_ _ ^ Y{A-ys )ys - {y.Y^ - 2{A + - y,)^)e cos 5 



(-y2 + (A-y,)2) 



3/2 



+ 



(r2 _ 2y2) (r - 2y,e cos (5)) iA^e cos 



y3 



(A3) 



r 



After solving the frequency response equation for the y- can be evaluated, thus determining whether or not para- 
resonance, Y and cos 5 can be found and the expressions metric excitation of the a;_-mode will occur. 
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